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Two non-equidistant grid implementations of infinite range exterior complex 
scaling are introduced that allow for perfect absorption in the time dependent 
Schrodinger equation. Finite element discrete variables grid discretizations provide 
as efficient absorption as the corresponding finite elements basis set discretizations. 
This finding is at variance with results reported in literature [L. Tao et al., Phys. 
Rev. A48, 063419 (2009)]. For finite differences, a new class of generalized Q-point 
schemes for non-equidistant grids is derived. Convergence of absorption is exponen¬ 
tial ~ Ax^~^ and numerically robust. Local relative errors < 10“® are achieved in 
a standard problem of strong-held ionization. 


I. INTRODUCTION 

In the numerical solution of partial differential equations (PDEs) for physical problems 
that involve scattering or dissociation one usually tries to restrict the actual computation 
to a small inner domain and dispose of the parts of the solution that propagate to large 
distances. The art of achieving this without corrupting the solution in the domain of interest 
is called to “impose absorbing boundary conditions”. Even without considering questions 
arising from discretization it is difficult to lay out a method that would provide perfect 
absorption in the mathematical sense. By perfect absorption we mean a transformation of 
the original PDE dehned through an operator D to a new one Da such that their respective 
solutions T and Tq agree in the inner domain and that Tq is exponentially damped in the 
outer domain. For reasons of computational efficiency we usually require Da to be local, i.e. 
composed of differential and multiplication operators, assuming that D is local. Without 
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this requirement one can often resort to spectral decomposition of D and apply the desired 
manipulations to each spectral component separately to obtain Da- However, this in general 
involves non-local operations with a large penalty in computational efficiency. 

Three local absorption methods for Schrodinger-like equations are particularly wide¬ 
spread in the physics community: complex absorbing potentials (CAPs), mask function 
absorption (MFA), and exterior complex scaling (ECS). In CAP one adds a complex po¬ 
tential, symbolically written as Da = D + Va, that is zero in the inner domain and causes 
exponential damping of the solution outside. The transition from the inner to the outer 
domain is smoothed to suppress reflections from the transition boundary. Such potentials 
are easy to implement and can be rather efficient, in particular if the spectral range that 
needs to be absorbed is limited. We restrict the dehnition of CAP to proper potentials, i.e. 
multiplicative operators that do not involve differentiations. CAPs of this kind are never 
perfect absorbers as defined above. 

MFA is arguably the most straight-forward idea: at each time step one multiplies the 
solution by some mask function that is smaller than 1 in the outer domain. In the limit 
of small reduction at frequent intervals this clearly approximates exponential damping in 
time. Also, as the solution propagates further into the absorbing layer, this translates to 
exponential damping in space. As in CAPs, the mask functions usually depart smoothly 
from the value of 1 at the boundary of the inner domain to smaller values, often to zero, at 
some hnite distance. In many situations MFA can be understood as a discrete version of a 
purely imaginary CAP 14 by dehning the mask function as 


Ma{x) = exp[-iAtVa{x)], 


( 1 ) 


where At is the time step. In such cases one finds similar numerical behavior for both meth¬ 
ods and the choice between MFA and CAP is only a matter of computational convenience. 

ECS is somewhat set apart from the first two methods in that it systematically derives Da 
from D by analytic continuation, trying to maintain the desireclproperties. If one succeeds, 
one obtains a perfect absorber in the mathematical sense l|-l3(]. This can be proven for 
stationary Schrodinger operators with free or Coulomb-like asymptotics and it has been 
demonstrated numerically for an important class of linear Schrodinger operators involving 
time dependent interactions j^. The method will be discussed in more detail below. 

In spite of being “perfect” and, as we will demonstrate below, also highly efficient, ECS 
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has remained the least popular of the three methods, although we may be observing a 
recent surge in its application jd-?]. The rare use may be related to the fact the ECS 
requires more care in the implementation than CAP and MFA. In fact, to the present date, 
the most efficient implementations have been by a particular choice of high order hnite 


elements (FEMs), named “inhnite range ECS” (irECS) j^, and other local basis sets such 
as B-splines js]. Such methods involve somewhat higher programming complexity than grid 
methods and, more importantly, pose greater challenges for scalable parallelization. 

When ECS is used in grid methods one usually introduces a smooth transition from the 
inner to the absorbing domain, sometimes abbreviated as smooth ECS (sECS, [9|). Reports 
about the efficiency of such an approaches appear to be mixed, usually poorer than in the 
FEM implementation, and certainly the results do not deserve the attribute “perfect”. We 
will state and demonstrate below what one should expect from a perfect absorber in compu¬ 
tational practice. The lack of perfect absorbers for grid methods is particularly deplorable, 
as grids are usually easily programmed and as they can also easily be applied to non-linear 
problems. 

In this paper we overcome the limitation of irECS to the FEM method and introduce 
grid-implementations that are comparable in absorption efficiency with the original irECS 
implementation, while maintaining the scalability of grid methods. We discuss two inde¬ 
pendent classes of grids: the FEM-discrete variable representations (FEM-DVR) and finite 
difference (ED) schemes. The FEM-DVR is a straight forward extension of FEM and con¬ 
trary to earlier reports in literature, it is fully compatible with ECS and irECS. 

For ED we introduce a new approach to non-equidistant grids which maintains the full 
consistency order of ED also for abrupt changes in the grid spacing. ECS on grids can be 
understood as a (smooth or abrupt) transition to a complex spaced grid in the absorbing 
region. The irECS scheme correspond to using an exponential complex grids for absorption. 

Apart from plausibly deriving the schemes, we demonstrate all claims numerically on a 
representative model of laser atom interaction. The computer code and example inputs have 
been made publicly available at [l0|. In particular we will show that the consistency order 
of the grid schemes is as expected ~ Ax^~^, where Q is the number of points involved in 
computing the derivative at a given point, and that Q can be increased to approach machine 
precision accurate results without any notable numerical instabilities. We will show that 
one can construct such schemes even when the solutions are discontinuous and that smooth 
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ECS bears no advantage over an abrupt transition from inner to absorbing domain. 

The paper is organized as follows: a brief summary of ECS is given and the FEM irECS 
implementation is laid out. Next we show that FEM-DVR is obtained from FEM simply 
by admitting minor integration errors. We will show that these errors do not compromise 
the absorption properties of irECS. The second part of the paper is devoted to the new 
ED schemes applicable for non-equidistant grids in general and for irECS absorption in 
particular. 


II. EXTERIOR COMPLEX SCALING 


Exterior complex scaling is a transformation from the original norm-conserving time de¬ 
pendent Schrodinger equation with a time dependent Hamiltonian H{t) and solution t) 
to an equation of the form 

r! 

( 2 ) 


where solutions Tq become exponentially damped outside some hnite region, while inside 
the hnite region the solution remains strictly unchanged: \l/(a:,t) = for |a;| < Rq 

and for all t. Although, to our knowledge, rigorous mathematical proof for this fact is still 
lacking, convincing numerical evidence for the important class of time dependent Hamil¬ 
tonians with minimal coupling to a dipole held has been provided [^. Apart from these 
fundamental mathematical properties, in practical application it is important that machine 
precision accuracy can be achieved with comparatively little numerical ehort. The partic¬ 
ular discretization scheme that provides this efficiency was dubbed “inhnite range exterior 
complex scaling” (irECS). In Ref. it is shown that, comparing with a popular class of 
complex absorbing potentials (CAPs), the irECS scheme provides up to 10 orders of mag¬ 
nitude better accuracy with only a fraction of the absorbing boundary size. This original 
formulation of irECS was given in terms of a hnite element discretization. To our knowledge 
(and surprise) irECS has not been implemented by other practitioners, although its good 
performance appears to have lead to a re-assessment of ECS methods for absorption and 
encouraging results were reported [#-[7|. 

In this section we give a brief formulation of ECS that will allow us to formulate the 
essential requirements for a numerical implementation. 
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A. Real scaling 


Exterior complex scaling derives from a unitary scaling transformation Ua from the orig¬ 
inal coordinates x to new coordinates y which is dehned as 

(f/.>!>)(!,) = A(!,)‘/>(A(!,)), Afe)- r \(s)ds (3) 

J —oo 

with a real scaling function 


0 < A(|/) 


1 for \y\ < Ro 
ag{y) for \y\ > Rq 


( 4 ) 


The transformation is unitary for any positive g{y) > 0 and any positive a. One sees that 
the transformation leaves {Ua'^){y) = invariant in the inner domain y < Rq and that 
it stretches or shrinks the coordinates for ag{y) ^ 1. Note that here X{y) only is required 
to be positive, but no continuity assumptions are made. 

Switching from the Schrodinger to the Heisenberg picture and considering Ua as a trans¬ 
formation of operators rather than wave functions, we can dehne a scaled Hamiltonian 
operator 

Ha = UaHUl (5) 


Clearly, as a unitary transformation Ua leaves all physical properties of the equation invari¬ 
ant. One can think of the transformation as the use of locally adapted units of length. 

Some caution has to be exercised, when \{y) is non-differentiable or discontinuous. Ob¬ 
viously, starting from a differentiable T the corresponding Ta will become non-differentiable 
or discontinuous to the same extent as \{y) is non-differentiable or discontinuous. As a 
result, we cannot apply the standard differential operators A or V to Tq. Of course, by con¬ 
struction we can apply the transformed = UaAU* and Vq = UaAU* to it. Conversely, 
the transformed —Aq and Va cannot be applied to the usual differentiable functions T, but 
only to functions obtained from differentiable functions by the transformation Tq = Ua"^- 
This simple observation will be the key to constructing numerically efficient discretization 
schemes for the scaled equations and also to write all transformed discretization operators 
in the simplest possible form. 
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A short calculation shows that the transformed hrst and 2nd derivatives have the form 


= i\{y) 


= -9{y) 
a 


hl>JJo 


= ^9(9)-‘'Vg(9)-‘V9fe)-‘'" 




( 6 ) 

( 7 ) 


here written in a manifestly Hermitian form. Potentials simply transform by substituting 
for the argument: 

Va{y) = V{A{y)). (8) 


B. Complex scaling 

Complex scaling consists in admitting complex a with A{a) > 0. To see how this leads to 
exponential damping of the solution, one can consider the transformation of outgoing waves 
/c > 0 at x —)■ cxo (assuming for simplicity g{y) = 1) 

T(x) ~ ^ '^aiy) ~ a^/2gifcy5R(a)-fcyO(a)_ 

Ingoing waves k < 0 would be exponentially growing and are excluded, if we admit only 
square-integrable functions in our calculations. This is the case if we calculate on a hnite 
simulation box [—L, L] with Dirichlet boundary conditions at ±L. 

For dehning 14 at complex values of a, there must be an analytic continuation of V{x) 
to complex arguments 14(|/) = V{A{y)) in the outer domain x > Rq. For being useful in 
scattering situations, the analytic continuation must maintain the asymptotic properties of 
the potential, such as whether it admits continuous or only strictly bound states. This is not 
guaranteed: for example, for arg(a) > 7r/4, complex scaling turns the harmonic potential 
from conhning oc x^ to repulsive cx —y‘^. No such accident happens in typical systems 
showing break-up or scattering with Coulomb or free asymptotics. A much more profound 
discussion of the mathematical conditions for complex scaling can be found in [n|. 

Apart from the exponential damping of the solutions, the second important property 
for application of ECS in time dependent problems is stability of the time evolution: the 
complex scaled Hamiltonian Ha{t) must not have any eigenvalues in the upper half of the 
complex plane. If such eigenvalues appear, they will invariably amplify any numerical noise 
as the solution proceeds forward in time and as a result the complex scaled solution will 
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diverge. Luckily, stability has been shown for a large class of Schrodinger-type equations. 
This includes time dependent Hamiltonians with velocity gauge coupling iV ■ A{t) to a 
time dependent dipole vector potential A{t). It is interesting to note that the length gauge 
brmulation of the same physical problem with the coupling f ■ S(t) in instable under ECS 
This may be surprising, considering that the two forms are related by a unitary gauge 
transformation. However, as the gauge transformation is space-dependent, its complex- 
scaled counterpart is not unitary and changes the spectral properties of the Hamiltonian. 

For a more detailed discussion as to why the solution remains invariant in the inner 
domain and under which conditions time-propagation of the complex scaled system is stable, 
we refer to [311 and references therein. 


III. THE FEM-DVR METHOD FOR IRECS 


Here we lay out how ECS is implemented for FEM-DVR methods. The FEM-DVR 
approach was introduced in 12j. Mathematically it differs from a standard hnite element 
method only by the admission of a small quadrature error. Therefore we hrst formulate 
the standard hnite element method in a suitable way. We limit the discussion to the one¬ 
dimensional case. Extensions to higher dimensions are straight forward and the problems 
arising are not specihc to the individual methods. 


A. A formulation of the finite element method 

In a one-dimensional hnite element method one approximates some solution T(a:), x G 
M piecewise on N intervals, the hnite elements [x„_i,a:„],n = With local basis 

functions f'^\x ),... /c = 1,..., that are zero outside the [a;„_i, a;„] one makes the ansatz 

N Q„ 

= ( 10 ) 

n=l k=l 

Note that interval boundaries Xn, number and type of functions can be chosen 
without any particular constraints, except for the usual requirements of diherentiability, 
linear independency, and completeness in the limit Qn oo. If the exact solution T 
is smoothly diherentiable, polynomials are a standard choice for However, at specihc 
locations or at large |a:| other choices can bear great numerical advantage, as will be discussed 




below. Equations of motion for the C are derived by a Galerkin criterion (in physics usually 
called the Dirac-Frenkel variational principle) with the result 




( 11 ) 


where the Hamiltonian H{c) and overlap S matrices are composed of the piece-wise matrices 


frin) _ 
^kl — 


qin) _ 
— 




W/ 


^n—1 




(n), 


( 12 ) 

( 13 ) 


Here and in the following we assume that the are real-valued. Complex functions used 
in practice, such as spherical harmonics, can be usually obtained from purely real functions 
by simple linear transformations. Clearly, we assume that H{t) is local, i.e. matrix elements 
of functions from different elements n ^ n' are = 0. Note that all basis sets that are related 
by a Qn X similarity transformation ^ire mathematically equivalent. 

These free parameters can be used to bring the to a computationally and numerically 
convenient form. 

So far, the ansatz admits ^(a:) that are not differentiable or even discontinuous at 
the element boundaries Xn- It is well known that for a correct dehnition of the discretization 
of the differential operators id^ and <9^ it is sufficient to ensure that the T {x) are continuous 
at Xn, if one secures that all operators are implemented in a manifestly symmetric form. The 


correct symmetric form is typica 


terms are dropped (see, e.g.. 
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ly obtained by a formal partial integration where boundary 
for a more detailed derivation). For example 




Oi 


(14) 


Note that explicit (anti-)symmetrization must also be observed for operators involving Erst 
derivatives, e.g. 


+ |)/,‘">) ^ t ((/i”’l94//">) - (94/f>|//”’)) . 


( 15 ) 


Continuity can be most conveniently realized by applying a similarity transformation 
on each element such that only the hrst and the last function on the element are non-zero 
at the interval boundaries and fixing these boundary values to 1: 



else = 0. 


( 16 ) 
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Implementation of these 2Qn conditions fixes only 2Q„ out of the Q\ free parameters in 
The remaining freedom can be used for further transforming the basis set. 

With such functions continuity can be imposed by simply setting equal the coefficients 
Cnk corresponding to the left and right functions at each element boundary Xn'- 

^n,Q — Cn+1,1 Vu. 

In the full matrices S and H{t) the corresponding rows and columns will be merged. One 
readily sees this amounts to adding and into S and H{t) such that the lower right 
corners of the n’th submatrix overlaps with the upper left corner of the n + 1st matrix (for 
an illustration, see, e.g. j^). 

In general, the will be full. One can use the remaining freedom in to bring 
the matrices to nearly diagonal form where only two non-zero off-diagonal elements 
remain and there are all I’s on the diagonal except for the hrst and the last diagonal entry. 
Complete diagonalization of S is inherently impossible without destroying the locality of the 
FEM basis. 

The non-diagonal form of S is the primary technical difference between grid methods and 
FEM. It is a signihcant drawback, in particular, when operating on parallel machines, where 
either iterative methods must be employed or all-to-all communication is required. This is 
not a problem of operations count; applying the inverse in its near-diagonal form with only 
two off-diagonal elements for each of the N elements can be reduced to solving a tri-diagonal 
linear system of size — 1. However, solving the tridiagonal system connects all elements 
to each other. In a parallel code where the elements are distributed over compute nodes this 
ensues costly all-to-all type of communication and may require complex coding, especially 
in higher dimensions. 


B. A formulation of the FEM-DVR method 

In FEM-DVR one reduces the overlap matrix to S' = 1 by admitting a small quadrature 
error in the computation of matrix elements. We introduce the FEM-DVR discretization 
based on the approach above. We choose our functions in the form 




( 18 ) 
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with polynomials of maximal degree Qn — 1 for pk{x) and a weight function v^'^\x). For 
such functions there is a Q„-point Lobatto quadrature rule 


Qn 




( 19 ) 


where the quadrature points include the interval boundary values x„_i = Si < Si < ... < 
^Qn = ^n- We can construct our basis functions using the Lagrange polynomials for the 
Lobatto quadrature points Sj as 


/F(-)=^n—. 

V^^){Sk) Sk - Si 


y(n)(^. 


X - Si 


( 20 ) 


which have the properties ffTTD . If instead of exact integration one contents oneself with 
(approximate) Lobatto quadrature one hnds a diagonal overlap matrix: 


Qn 

~ 5^WiPfc(si)Pi(si) = WkSki. (21) 

i=l 

In fact, exact integration is only missed by one polynomial degree, as Lobatto quadrature is 
exact up to degree 2Qn — 3, while our are degree Qn — 1- The two degrees lower accuracy 
of Lobatto compared to standard Gauss quadrature is the penalty for hxing si and sq^ to 
coincide with the interval boundaries. 

A further advantage of FEM-DVR over FEM is that Lobatto quadrature is applied for 
all multiplicative operators, not only the overlap. By that all multiplication operators are 
strictly diagonal and allow highly efficient application. The advantage is mostly played out 
in higher dimensions, where the exact basis set representation of a general potential would 
be a full matrix. Derivative operators are full in FEM-DVR, but usually they come in the 
form of a short sum of tensor products, which again can be implemented efficiently. 


C. ECS and irECS for a FEM-DVR grid 


The favorable absorption properties of ECS in general and of the particular implemen¬ 
tation by irECS were hrst reported in Ref. and have since be used to solve several 


challenging problems in the strong laser-matter interactions 14l-ll8|. All these calculations 
were performed using a FEM basis. 


In fact, in Ref. 19| severe instabilities were reported for ECS absorption using FEM- 
DVR discretization for a simple test-problem where the irECS showed perfect absorption. 
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f B 


In Ref. m we speculated that the the approximate quadratures inherent to FEM-DVR were 
to blame. This is a plausible possibility, as analyticity plays a crucial role for ECS and small 
integration errors by using Lobatto quadrature instead of evaluating integrals exactly might 
destroy perfect absorption. Now we show that this speculation was incorrect, that FEM- 
DVR gives numerical results of the same quality as FEM, and that the problems encountered 


m 


19| must have had a different origin. 


All calculations below were performed using the tRecX code, which together with the 
relevant example inputs has been made publicly available 101. 


1. Model system 


In all numerical examples in this paper we use as a model Hamiltonian the “one¬ 
dimensional Hydrogen atom” in a laser held (using atomic units h = e = rrie = 1) 

H(t) = - ,A(t)d, - (22) 

^ + 2 


where A{t) is the laser held’s vector potential (in dipole approximation). Remarkably, the 
ground state energy of the system is exactly at -1/2, as in the three-dimensional Hydrogen 
atom. It has been demonstrated that the mathematical behavior of absorption in this simple 
system generalizes to the analo gou s Schrodinger equations for one- and two-electrons systems 


in up to 6 spatial dimensions 
form 
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, For all studies below use a vector potential of the 

(23) 


TTt 


A{t) = Aocos^(^) sin(a;t) for t G [-T,T] 


with Aq = 1.3 aw, uj = 0.057, and T = Ott/o;. In more physical terms, the parameters 
translate into a pulse with central wavelength of 800 um, peak intensity 2 x lO^^hF/cm^, 
and a FWHM duration of three optical cycles. Such a pulse depletes the initial ground 
state of the system by about 50%, which is all absorbed at the boundaries. In this type 
of strong-held ionization processes emission occurs over a very wide spectral range. At our 
parameters outgoing wave-vectors cover the whole range from zero up to ~ 2 au before 
amplitudes drop to below physically relevant levels. This broad range of outgoing wave 
vectors poses a particular challenge for absorption. At more narrowly dehned ranges of 
wave vectors absorption can be achieved by a variety of methods by tuning the parameters 
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to the specific wave vectors. One of the advantage of ECS is that it can be applied over the 
whole range without the need to adjust to the particular form of emission. 

An important physical parameter of strong field photo-emission is the “quiver radius”: a 
classical free electron will oscillate with an amplitude oq = Aq/o; in the laser field. At our 
parameters one computes a quiver radius of oq ~ 23. This gives a rough measure for the 
radius up to where one needs to preserve the solution without absorbing and it motivates 
our choice of Rq > 25. Note that, if so desired, ECS allows choosing arbitrarily small Rq 
(including Rq = 0) such that flux may propagate deeply into the complex scaled domain and 
return to the inner domain without necessarily corrupting the solution (see Ref. j3| for more 
details). This fact further corroborates that, mathematically speaking, ECS is a lossless 
transformation. The loss of information by the exponential damping is purely numerical 
due to the limited accuracy of any finite representation of the solution. 


2. Implementation of ECS 


We use the simplest scaling function g{y) = 1, i.e. 


X{y) = a for \y\ > Rq. 


(24) 


The scaled Hamiltonian is 


Ha{t) 


-^[dl]a - iA{t)[dy]a 


1 


with the scaled derivatives 


[dy]a, [0% 


on\y\<Ro 
^dy, on \y\ > Rq. 


(25) 


(26) 


The scaled solution will have the form UaAl for some differentiable T. In particular, 4/^ has 
a discontinuity by the factor when crossing the scaling radius Rq. In a finite element 
scheme it is easy to implement such a discontinuity: we choose two element boundaries to 
coincide with the lower and upper boundaries of the inner domain RRq = yn±- Then all 
functions on the outer domain are multiplied by 




for yn-i > or < -Rq. 


(27) 
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The desired discontinuity is ensured by equating the coefficients corresponding to the bound¬ 
aries just as continuity is ensured at all other boundaries. Conditions on the derivatives 
can be omitted for the same reasons and with the same precaution about using explicitly 
symmetric forms of the operators as discussed above, Eqs. d and (1T51) . 

In practical implementation, multiplying the function translates into a multiplication of 
the element matrix blocks and by a: 




e{n) 






for y„ < -Ri> or flo < Hn-I- 


(28) 


Obviously, the only non-trivial effect of this extra multiplication by a appears at blocks to 
either side of ±i?o- Also note that, while the overlap matrix blocks remain unchanged 
except for the multiplication by a, one must use the properly scaled operators Ha{t) for 
evaluating the scaled matrix blocks 

Complex scaling now means that the operator is analytically continued w.r.t. a. There 
is a seeming ambiguity as how to deal with complex conjugation of in the scalar 

products. One might suspect that in fact there should be |a| appearing as a factor for 
the matrices rather than a. Clearly, this would pose a problem as the modulus is not 
an analytic function and analytic continuation of the operator would be doomed. Closer 

Tom the ket functions 


inspection shows that bra functions 
l"^), exactly such that the conjugat: 
appearing in the operator that are extended to complex values 


I must be chosen differently 
, exactly such that the conjugation of jg undone, see j^, 3, 


2l| . Thus, it is the a 


3. irECS discretization 

The irECS version of ECS greatly enhances computational efficiency by replacing the 
Dirichlet conditions at the finite boundaries ±L with a computation on an infinite interval 
where exponentially decaying basis functions ensure decay —)> 0 as ||/| —?■ cx). For our two- 
sided infinity this amounts to formally choosing —|/o = = oo and using the weight 

functions v{y) = exp{±'jy) at the first and last interval, respectively. The finite inner 
domain of the axis is divided into elements of equal size. We construct the as in (lT 6 l) . It 
has been investigated earlier how errors of irECS in the FEM implementation behave with 
order, number of elements, complex scaling angle 6, scaling radius, and exponential factor 
exp(— 7 I 1 /I) on the infinite intervals. Summing up those results, irECS absorption is highly 
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efficient and, within reasonable limits, quite insensitive to these details of the discretization. 
For a quali- and quantihcation of this statement we refer to Ref. j^. 

The irECS idea readily carries over to FEM-DVR, if we use a Radau quadrature formula 
for the inhnite intervals. Radau quadrature hxes only one quadrature point at the hnite 
left or right boundary of the interval. One adjusts the remaining quadrature points and the 
quadrature weights for the specihc weight function |u(|/)P such that with quadrature 
points integrals become exact up to polynomial degree — 2. 

For demonstrating that the numerical behavior of FEM-DVR and FEM are equal for all 
practical purposes, we use a hxed set of discretization parameters. We choose the complex 
scaling parameter a = exp(i/2) and radius Rq = 25 with 10 intervals of size 5a.u. in the 
inner domain [—Ro,Rq\. For the inhnite element basis we use the v{y) = exp(—1|/|/2), i.e. 
7 = 0.5. We use a uniform order Qn = Q in the inner domain for all elements and two 
inhnite elements of order Qi = Qj^f =: A in the outer domain to either side. We study the 
variation of the results with Q and A. The same functions with Lagrange polynomials 
at the Lobatto points (hnite elements) or Radau points (inhnite elements) are used in FEM 
and FEM-DVR. For FEM-DVR, this choice is by dehnition. For FEM the exact choice of 
the polynomials is unimportant: results near machine precision can be obtained with any 
set of polynomials, if only one avoids ill-conditioning problems as they typically arise in too 
simplistic choices, such as monomials. In fact, in all previous calculations we had derived 
our basis from Legendre (hnite) and Laguerre (inhnite range) polynomials. For FEM, we 
compute all matrix elements to machine precision using a recursive algorithm. In the given 
basis, FEM-DVR simply consists in replacing the exact integrals with Q-point Lobatto and 
R-point Radau quadratures on the respective elements. 

Throughout this work we assess the accuracy of the solutions by computing the local and 
maximal relative errors of the probability density at the end of the pulse p{x) := |\l/(x, T)p: 


e(x) 


p{x) - pojx) 
Po{x) 


Co := max e(x). 
xG[— i?o,i?o] 


(29) 


The reference density po is drawn from a large, fully converged calculation. 

For time-integration we use the classical Runge-Kutta scheme with step-size control based 
on the maximal error of the coefficient vector components Cnk- This universally applicable 
method was selected to facilitate comparisons between the methods, without any attempt 
to optimize its performance. 
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FIG. 1: Density p{x) (upper panel) and discretization error e(x) at the end of a 3-cycle laser pulse 
at 800 nm central wave length and peak intensity 2 x 10^^. The e(x) are barely distinguishable, 
solid line: FEM, dashed: FEM-DVR. The more accurate result at e(x) < lO”'^ was obtained with 
discretization {Q,A) = (12,18), 186 points, errors 10 ^ are reached with {Q,A) = (9,15), 150 
points. (See text for the definition of Q,A and e(x).) 

Figure [T] shows po{x) and the e{x) corresponding to two pairs of FEM and FEM-DVR 
calculations with errors e{x) < 10“^ and e{x) ^ 10“^. One observes that FEM and FEM- 
DVR produce, at equal discretization size, equally accurate absorption with no obvious 
accuracy advantage for either method. 

Fig. [2] shows thee discretization errors in the inner domain and complex scaled domain 
independently. Computations are with Q, Aq varying Q for the inner domain and Qq, A 
varying A for absorption. The joint parameters Qo, Aq are chosen such that the solution is 
maximally accurate. As expected, errors drop exponentially with Q, and also absorption 
improves exponentially beyond a minimum number of A > 12. 
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FIG. 2; Convergence of the peak error eo- Left: as function of discretization order Q in the inner 
domain keeping A = Aq fixed, right; as a function of absorption order A, keeping Q = Qq fixed. 
Solid lines: FEM, dashed FEM-DVR. System parameters as in[TJ 

D. Discussion and conclusions on FEM-DVR 


The simple conclusion of this first part is that FEM-DVR is completely at par with FEM, 
at least as far as perfect absorption is concerned. Considering the great simplihcation for 
the computation of integrals and, more importantly, the gains in operations count, ease of 
implementation, and ease of parallelization it is certainly to be preferred over a full blown 
FEM method. There is only a single point where we see some advantage of computing 
the integrals exactly: FEM eigenvalue estimates are variational upper bounds, while FEM- 
DVR approximations may drop below the true value. The actual errors are, according to 
our observation, of similar magnitude in both methods. In practice, the upper bounding 
property will rarely be of great importance. 

Our conclusion is at variance with previous reports on absorption using FEM-DVR. In 


Ref. [1^ the inhnite range idea was not used. However, at some expense as to efficiency, 
also a standard DVR discretization using only polynomials without exponential damping 


produces highly accurate results. This was reported for FEM in and can be reconhrmed 
for FEM-DVR here. Fig. [3] shows the errors e{x) of a discretization where the inhnite 
intervals were replaced by an increasing number of hnite elements of the constant order Qq. 
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FIG. 3; Local error e{x) for absorption in a finite box as the box [—Xmaxj^max] is extended by 
adding elements in the absorbing domain. Size of all elements is Xn — Xn-i = 5 and order is fixed 
at Q = ^4 = 21. 


One sees that errors can be very well controlled and no artifacts 
I a; I ~ i?o- In spite of similar discretization sizes, the calculations in 
instable, especially near |a:| = Rq. 


appear also in the area 
19| were reported to be 


IV. GENERALIZED ED SCHEMES 

In FD schemes there is usually no explicit reference made to an underlying basis, but 
rather one is contented with representing the wavefunction at the grid points. In reality also 
FD uses a hypothesis for evaluating the derivatives: in standard applications one assumes 
that in the vicinity of a grid point Xj the solution can be well approximated by a linear 
combination of polynomials fl^\x). We hrst re-derive the standard schemes in analogy 
to the discussion of the FEM method above and then generalize the approach for non- 
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differentiable solutions. 

Standard symmetric FD schemes on equidistant one-dimensional grids Xj = Xq + jh,j = 
0,..., J are obtained by assuming that in the vicinity of each point Xj one can write 

p 

~ fk\x)ck- (30) 

k=-P 


The “interpolation hypothesis” is that the functions are polynomials. We treat here 
only the case of symmetric schemes with an odd number of Q = 2P + 1 points and use the 
notation := f^^\xj+i) with the index ranges /, /c G {—P,... ,P}. The coefficients Ck can 
be obtained from the neighboring function values := "^{xj+i) as 


p 

Ck=^ 

i=-p 


pU) 


-1 




j+i- 


kl 


One readily obtains an approximation to the derivative {dx'^){xj) =: T'- as 


(31) 


~ 'PiSjiPxj) 

l,k 


pU) 


=: 


j+l 

kl i^_p 




j pi¬ 


rn) 


Arranging the pointwise hnite difference rules into a matrix = dP we hnd for the 

hnite difference approximation of the hrst derivative 




(33) 


The same construction principle can be used for higher derivatives or actually any operator 
composed of derivatives and multiplicative operators. 

Eq. fl3^ is suitable for the construction of the schemes in numerical practice, if only one 
avoids ill-conditioning of P^^K Almost any choice for the fP, e.g. orthogonal polynomials 
or Lagrange polynomials at well-separated support points, will suffice. 


A. Non-equidistant, discontinuous, and complex scaled schemes 

A standard way of constructing FD schemes for non-equidistant grids is based on the 
very same coordinate scaling discussed in the preceding section. Let us assume that the 
non-equidistant grid is dehned by some monotonically increasing function as 


Xj = A{yj) with yj = 0, h, 2h ,..., 


(34) 
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which transforms the representation on the non-equidistant Xj grid into a representation on 
the equidistant Hj grid. For example, for exponential sampling of a; > 0 one can choose 
A(|/) = exp(?/) — l,j = 0,..., iV — 1 oi y = log(l + |a;|). For deriving the necessary 
transformation of the operators and for constructing norm-conserving schemes, it is useful 
to adhere to the notation of section and consider the coordinate transformation as a 
unitary transformation d' —d'a = tAd', Eq. ([3]). 

The key to suitable FD schemes in the transformed coordinates is to realize that the 
interpolation hypothesis for the transformed solution in general must differ from the original 
one. Suppose a set of functions is in some sense an optimal interpolation hypothesis for 
the unstretched T(a;) near the point Xj. Then the transformed functions 

/“(») - (35) 

will be equally optimal for interpolating Ta(|/) near the point yj = A~^{xj). 

Assume that the interpolation hypothesis (x) of the original problem are polynomials, 
i.e. the standard finite difference scheme of a given order. Then in almost all cases the equally 
optimal interpolation hypothesis f^liy) = X{y)~^^'^f^^\A{y)) for the transformed solution 
Ta(?/) will not he polynomials and one needs to re-derive the corresponding finite difference 
scheme by Eq. fl32l) . 

In practice, the procedure for obtaining schemes for any transformed linear operator 
Ba = UaBU* is very simple. Analogous to Eqs. fl30|) and fl3^ one writes 

p p 

'SM = (t?.!')(!/)« ^ {ujh(y)ct = HyY'VX’(Hy))ct (36) 

k=-P k=-P 


P 

(B.'ta)(9j) = (!/„BC/:)4'.(9 j) « (UaBfY>)(y)ck 

k=-P 

P 

= E Hyiy'YBYA{A{yj)W 

k=-P 


( 37 ) 


With the notation 



= \{yiY/‘^{f]^\A{yi)) one obtains the adjusted scheme at point 


Vj as 



Q 

EA(ft)‘'"(B/f)(A(!/j)) 


k=l 



( 38 ) 
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For the first derivative [dy]a we insert = 


Aj) 

Jk 


into Eq. fl38|) . for the second derivative 


[dy]a we insert etc. As long as all transformations Ua are unitary, one obtains 

symmetric schemes for \l/a, provided that the original problem gives symmetric schemes. 

The adjustment of is particularly important, when the transformation A is not ana¬ 
lytical. This is for example the case, when we want to switch from constant spacing Axi in 
one region to a different constant spacing Ax 2 in some other region, possibly with a smooth, 
say linear, transition in between. If we attempt to approximate the transformed solution 
Ta(?/) by higher order polynomials in y, i.e. when we use standard higher order FD schemes, 
we will observe a loss of the convergence order. Assuming the original T (x) has a convergent 
Taylor expansion, the solution w.r.t. to the transformed coordinate A^/^(|/)Ta(|/) = \1/(A(|/)) 
only has a continuous hrst derivative dy. Already the second derivative will be discontinu¬ 
ous and all higher derivatives show (j-function like singularities at the boundaries between 
constant and linearly changing spacing, which disqualihes any attempt of expanding into a 
convergent Taylor series. Making the transition smoother only postpones the problem to 
higher orders. On the other hand, we will demonstrate below that adjusting and using 
Eq. fIMll for the operators preserves the approximation order and no extra computational 
cost ensues apart from the initial construction of the local schemes 

From the discussion it is clear that we can also use transformations generated by non- 
differentiable A(x), as they arise in ECS. All we need to do is to replace polynomials of 
standard FD schemes with their transformed counterparts. Finally, also the idea underlying 
the irECS discretization discussed above can be transferred to FD schemes. Considering the 
approximately exponentially spaced Radau quadrature points appearing in the FEM-DVR 
implementation of irECS, we see that the essential point of irECS is that the function is 
sampled at rapidly increasing spacing rather than uniformly. By uniting complex scaling 
with non-equidistant sampling we will obtain nearly as efficient absorption scheme with FD 
as irECS with FEM-DVR discretization. 


cU) 


B. Absorption with equidistant FD grids 


We hrst investigate complex scaling with equidistant grids for the model Hamiltonian 
fl22ll . We will show that comparable accuracies are reached for FEM-DVR and FD at equal 
grid sizes and equal orders 
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We again use the transformation defined by Eq. fl24|) with a = exp(i/2) and Rq = 25 and 
FD schemes are constructed according to Eq. fl38|) . The same constant order = 21 is used 
in FD and FEM-DVR calculations on 1200 grid points in a computational box [—150,150]. 
Fig. m compares complex scaled densities |T(j(a:)p and errors of the two methods relative 
to a converged calculation on a large box without using absorbing boundaries. Comparing 
unsealed with complex scaled |\['a(x)| densities one can clearly discern the suppression 

of the density to below 10“^^ as the solution approaches the box boundaries. On the level of 
densities FD and FEM-DVR results are indistinguishable. Also, relative to a fully converged 
unsealed calculation, errors in the inner region are comparable. 

For completeness we show in Fig. [5] that the error drops exponentially with the FD order 
Q. This proves the full consistency and that the non-differentiable nature of the complex 



X (a.u.) 


FIG. 4: Densities of the unsealed |^(x)p (dot-dashed line) and the complex scaled |'I'a(a^)P system 
(solid line). Damping by complex scaling at large |x| is clearly exposed. Inset: relative errors e(x) 
in the inner domain for DVR and FD at equal grid sizes and order. 
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FIG. 5: Convergence of peak relative error eo on the inner domain: (a) convergence with Q for 
1600 equidistant grid points on [— 200 , 200 ], (b) convergence with Ax, dashed lines Ax^“^. At 
Q = 15 the actual convergence is only ~ Ax^^. 

scaled solution is fully accounted for by the generalized FD scheme. 


C. Absorption with exponential FD grids 


Having established standard ECS for FD grids, as a last step we implement the idea 
underlying irECS also for FD schemes. As discussed, this consists in exponentially sampling 
the scaled solution \l/a in the absorbing domain. In the notation introduced above, exponen¬ 
tial sampling can be achieved by applying a unitary transformation to the complex scaled 
solution 

= (39) 


where Wj is a real scaling transformation with 


A(2;) 


1 for \z\ < Rq 

exp[7(|2;| — i?o)] for kl ^ Ro 


( 40 ) 


For the comparisons we use a hxed complex scaling radius Rq = 25 and hxed complex 
scaling phase a = exp(i/ 2 ) as above. All complex scaled calculations are performed on the 
interval [—100,100] with a total of 400 grid points. This allows accuracies of the density of 
< 10“^ in unsealed domain [—25, 25], if absorption is perfect. Fig.| 6 ]shows the complex scaled 
|'ha(|/)p with equidistant grid and the complex scaled [T-y^a(.s)[ with exponential grids in the 
scaled domain for 7 = 0.1,0.2,0.35. The spatial damping of the solution by the complex 
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scaling transformation is clearly visible. Near the boundaries the equidistant grid solution 
|\l^aP drops to ~ 10“®. The effect of the Dirichlet boundary condition at \y\ = 100 is clearly 
visible, but has no consequences on the accuracy level of interest. The exponential grid 
maps the absorbing domain into smaller boundary layers outside [—25,25]. An optimum is 
reached near 7 = 0.2: the density drops to below the level < 10“® for all \z\ > 40. Larger 
contraction to 7 = 0.35 does not lead to further gains: although the solution shrinks to 
smaller \z\ initially, the exponential spacing Xj = x{zj) is becoming too wide to represent 
the solution in the asymptotic region. An artefact on the level ~ 10“^ appears which causes 
reflection errors. As shown in the hgure, the artefact can be fully suppressed by reducing 
the grid spacing, but this defeats the purpose of minimizing the number of points used for 
absorption. 

Finally, we investigate the dependence of absorption on 7 and 6 for FD and FEM-DVR. 
In both methods we use the same number of 201 discretization points in the unsealed domain 
[—25, 25] and the same order Q = 21, i.e. a 21 point scheme for FD and degree 20 polynomials 
for FEM-DVR, and we use the same number of A = 60 points for absorption. At the grid 



FIG. 6 : The irECS method for FD schemes. Solid lines: equidistant grid 7 = 0 and exponential 
scaling in the absorbing region. Spacing of the transformed grid Az = 0.25 and order Q = 21 . At 
7 = 0.35 insufficient sampling leads to artefacts. Dashed line: with spacing Az = 0.125 exponent 
7 = 0.125 does not show artefacts. 
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spacing of Az = 0.25 this amount to Dirichlet boundary conditions at \z\ = 40 for FD. 
Note that this results in matrices the same size and of comparable sparsity in the inner 
region, with band-width 21 for FD and near block-diagonal matrices of blocks 21 x 21 for 
FEM-DVR. Fig. [7| shows the relative errors of |4/(a:)p in the inner domain as obtained with 
either method. In FD the range of admissible 7 remains smaller than in FEM-DVR: as large 
7 lead to a stronger contraction of the wave function, we see that the profit from the irECS 
idea is somewhat lower in FD. 


D. Smooth exterior scaling 


We have demonstrated that the generalized FD schemes as well as FEM-DVR grids allow 
for discontinuous scaling functions A(a;). A fortiori we expect smooth exterior scaling to 
work correctly, if implemented by the above principles. This is indeed the case. We use the 
scaling function 


Ky) = < 


1 

1 + (n 


a 


for \y\ < Rq 

l)s(|l/|) for Ro < \y\ < Ri 
for \y\ > Ri, 


(41) 


\ 


where the 3’’'^ order polynomial s{y) smoothly connects the inner domain with the region 
\y\ > Ri such that X{y) is differentiable at y = Rq and y = Ri. Applying the same procedure 




FIG. 7: Relative error of FEM-DVR (left) and FD (right) schemes as a function of exponent 7 
and scaling angle 6 at equal number of discretization points and comparable matrix sparsity. The 
admissible 7 -range is narrower for FD, 9 can be chosen in a wide range with either method. 
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as above and using Ri = 45, i.e. smoothing over a range of 20 au, we find essentially identical 
results as for the abrupt transition used earlier, see Fig. |8l 

The usual motivation for smooth scaling is not to use the generalized schemes but to 
apply standard schemes. There is no unique way for defining such a scheme for a polynomial 
interpolation hypothesis. There are two different possibilities: one by bringing the scaled 
second derivative Eq. ([7]) to a form that allows the use of standard FD schemes for dy and 
dy without increase of band-width 



Uadlu: 


2X^^^X^ X^^2X^ 4 VaV 


(42) 


In a second approach, one can use the procedure for constructing finite difference schemes 
that lead to Eq. fl32ll for the complete operator fl42ll with polynomial interpolation in place of 
the correctly scaled polynomials. The latter fully uses the polynomial interpolation hypoth- 



X (a.u.) 

FIG. 8: Abrupt (solid line) vs. smooth ECS (dashed line). Exponential snppression is slightly less 
with smooth scaling, errors in the unsealed region are identical. 
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esis, the former introduces some additional approximations in evaluating the transformed 
derivatives of the polynomials [d‘^aPk{y)- Note that in both cases the polynomial interpola¬ 
tion hypothesis is unjustihed and the approximation cannot be successful beyond the lowest 
orders. 

For the present demonstration we only study the second possibility, i.e. using the full, but 
incorrect polynomial interpolation without making further approximations. The limitation 
of convergence order for the polynomial interpolation hypothesis applied to a smoothly 
exterior scaled grid is illustrated in Fig. |9l For orders Q=3 and Q=5 polynomial and 
scaled interpolation give essentially the same result. However, while the correctly scaled 
interpolations shows steady exponential convergence at Q = 7, 9, no further gains can be 
made by increasing the order with the unsealed polynomial hypothesis. For higher accuracies 
one is forced to increase the of number grid points, which causes exceeding numerical cost 
not only in terms of the problem size but also in stiffness of the equations. 

Without demonstration, we only remark that also for the FEM-DVR scheme smooth 
transition does not bear any advantage over the abrupt transition used in the FEM-DVR 
calculations above. On the contrary, the rather complicated smoothly scaled operators must 
be programmed and the hnite modulation of the solution in the transition region usually 
requires more grid points there, which in turn may raise the stiffness of the dynamical 
equations. 


V. CONCLUSIONS 

With the present studies we have demonstrated that grid methods allow for highly efficient 
absorption schemes. The irECS method had been originally developed using a hnite element 
implementation, where also the far superior performance of irECS compared to multiplicative 
CAPs and MEA had been highlighted. With the present extension of the method to FEM- 
DVR and FD grids one has discretizations that are easy to implement, have low hoating 
counts, and are straight forward for parallelization. The latter point may be the most 
important advantage. While in transferring irECS to FEM-DVR no practical problems of 
any kind were encountered, we had to derive a new approach to FD for complex scaled 
problems. In fact, the computation of scaled schemes is technically simple and provides 
schemes that are as nearly as efficient in terms of discretization size as the FEM and FEM- 
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FIG. 9: Convergence with order Q for scaled and standard polynomial schemes. Smooth exterior 
scaling with a third degree polynomial smoothing function on an equidistant grid were used. Rel¬ 
ative errors e{x) in the unsealed region [—25, 25] are show for Q = 3,5, 7,9. Scaled schemes errors 
decrease exponentially with Q (dashed lines). For polynomial schemes (solid lines) no gain is made 
for Q > 7 ). Calculations are with 600 points on the interval [—80,80]. The thick solid line is for 
a standard polynomial scheme with Q = 7 and as many as 4800 points on [—80,80]. 

DVR schemes of the same order. In particular, they can be pushed to high orders with 
exponential reduction of errors to near machine precision. 

We have shown that irECS can be considered as ECS with an exponential grid in the 
absorbing domain. Both, in FEM-DVR and in ED the transition from the equidistant 
unsealed to the exponentially spaced absorbing region is numerically seamless, i.e. it does 
not produce any artefacts compared to an equidistant grid. The numerical gain by the 
reduction of grid points is substantial and no increase in stiffness was observed. 

We have further studied smooth exterior scaling and shown that one can smoothen the 
















transition between scaled and unsealed region. This case was only discussed for the FD 
implementation, as one usually resorts to such a procedure to circumvent manifest problem 
that arises for using standard FD schemes across abrupt transitions. Absorption works 
for smooth scaling as well as it does for abrupt scaling, but only if scaled, non-polynomial 
schemes are used. With standard polynomial schemes, lack of differentiability leads to a 
strict limitation of the consistency to the point where the corresponding higher derivatives 
of the smoothing function become singular. Smoothing does not bear any mathematical or 
computational advantage as compared to abrupt transitions, on the contrary, it tends to 
complicate implementation as the essentially arbitrary smoothing transformation must be 
incorporated into the scheme. Therefore, at least from the present work, we would advice 
to use abrupt changeovers wherever possible. 

The FD schemes for non-equidistant grids, in particular the abrupt and reflectionless 
transition between grid spacings occurring in abrupt exterior scaling, may be of broad inter¬ 
est for uses beyond the problem of absorbing boundaries considered here. Locally adapted 
FD grids are frequently used in literature. To the best of our knowledge, the approach 
presented here and proven to fully maintain convergence orders is novel. In future work, 
we plan to investigate the problem for non-equidistant grid representations of Maxwell’s 
equations. 
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